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ABSTRACT 

We develop a new code, the Hierarchical Bound- Tracing (HBT for short) code, to 
find and trace dark matter subhaloes in simulations based on the merger hierarchy 
of dark matter haloes. Application of this code to a recent benchmark test of finding 
subhaloes demonstrates that HBT stands as one of the best codes to trace the evo- 
lutionary history of subhaloes. The success of the code lies in its careful treatment of 
the complex physical processes associated with the evolution of subhaloes and in its 
robust unbinding algorithm with an adaptive source subhalo management. We keep 
a full record of the merger hierarchy of haloes and subhaloes, and allow growth of 
satellite subhaloes through accretion from its "satellite-of-satellites" , hence allowing 
mergers among satellites. Local accretion of background mass is omitted, while re- 
binding of stripped mass is allowed. The justification of these treatments is provided 
by case studies of the lives of individual subhaloes and by the success in finding the 
complete subhalo catalogue. We compare our result to other popular subhalo finders 
and show that HBT is able to well resolve subhaloes in high density environment and 
keep strict physical track of subhaloes' merger history. This code is fully parallelized 
and freely available upon request to the authors. 
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1 INTRODUCTION 

In the hierarchical universe, cold dark matter haloes grow 
mainly through mergers with surrounding smaller haloes. 
After the merger, the imprints of progenitor haloes are not 
wiped out but in fact they can survive for quit e a long time 
as se l f-bound substructu r es called subhaloe s dMoore et al 
19981: iGhigna et all 1 19981 : iKlvpin et al] Il999l : 



Moore et al 



19991 ). Galaxies form inside dark matter haloes and can be 



traced by dark matter subhaloes after the merger. Because 
the non-linear growth of structure in the dark matter com- 
ponent can be well produced in N-body simulations, it has 
become a standard approach to build galaxy formation mod- 
els on top of the dark ma tter halo merger history (see e.g, 
lBaughll2006l : iBensonllioiol . for recent reviews). Constructing 
the full hierarchy of the merger history requires the identi- 
fication and linking of dark matter subhaloes across cosmic 
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time. Many algorithms have been developed to accomplish 
this job, all based on some of the following characteristics of 
a subhalo: 

(i) It is an overdense region inside its host halo 

(ii) It is self-bound so that it is dynamically significant 

(iii) It was a halo before it merges into its current host 
halo 

Most subhalo finders utilize only contemporary parti- 
cle distribution and focus on the first two characteristics, 
while the merger tree is constructed by a subsequent match- 
ing of subhaloes at different epochs . For example, the Hier- 
archical Friends-of-Friends (HFoF: [Kh^pin et al.lfl999l )) al- 
gorithm, which is an extension t o the standard Friends- 
of-Friends (FoF: iDavis et allllQSSl ) halo finder with multi- 



nly . 

AMIGA Ha lo Finder (AHF: iKnoUmann fc Knebd liooi ). 
SUBFIND (ISpringel. Yoshida. fc Wh"it3 1200 ll ) and SKID 
iGhigna et al. 19981 ). which collect local overdense parti- 
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cles and then eliminate unbound particles, are based on the 
first two characteristics. Although these percolation based 
algorithms are able to find subhaloes using a single simu- 
lation output, a strong resolution problem can arise in the 
central part of the host halo due to ambiguity of separat- 
ing member particles of a subhalo from the background 
ones in the high density region jGill. Knebe. fc GibsonI 
12004 iMuldrew. Pearce. fc Poweij I2OIII '). In fact in a re- 
cent extensive halo-finder comparison project (|Knebe et alj 
I2OIII I. those participating finders based on configuration 
information only all struggle to recover substructures in 
the central high-density region of a host halo (see also 
13. f[) . One way to improve the resolution of subhaloes in 
the high density region is to use six dimensional phase 
space information, as done in, e.g, the Hier archical Struc- 



ture Fi nder (HSF: iMacieiewski et all [20091 )') and ROCK- 
STAR (|Behroozi et al.ll201ll ). Another way out of this prob- 
lem is to appeal to the third characteristic and utilize 
the evolutionary history of subhaloes. Because subhaloes 
are remnants of dark matter haloes, they can be identi- 
fied by tracing the member particles of their pr ogenitor 
haloes . The first attempt along this direction was iTormenI 
lll997^ who just considered the th ird characteristic. Later 
in lTormen. Diaferio. fc Svej (|l998l ) self-boundness was also 
added to define a subhalo. Only those haloes which fall di- 
rectly i nto the final halo of i nterest were examined until the 
work of lGiocoli et al.l (|201[ ) (hereafter GfO) where they ex- 
tended their SURV code ( Tormen. M oscardi ni. fc Yoshidal 
I2OO4I : iGiocoh. Tormen. fc van den Boschi i200ij ') to include 
subhaloes inside subhaloes. The same method was also 
implemented in the 'ML APM halo tracker' (MHT; 
iGill. Knebe. fc Gibson|[20o3 ). 

While simple as the idea looks, the subhalo identifi- 
cation through tracing a merger history still faces several 
difficulties in practice. The most challenging aspect of the 
problem is how to trace subhaloes robustly over several or- 
bital periods. Mass loss is the primary process associated 
with the subhalo evolution due to gravitational striping and 
harassment. In such a simplified picture, once a particle be- 
comes unbound to a subhalo, it no longer needs to be traced 
any more. However, as we will show, this kind of successive 
tracing is dangerous since artificial loss of bound particles 
can accumulate through every tracing step, eventually trig- 
gering a runaway loss of a subhalo's bound particles. On 
the other hand, if one always traces all the particles from a 
subhalo's progenitor halo, a straight-forward unbinding al- 
gorithm may also fail to find a self-bound structure once the 
subhalo has been substantially stripped compared to its pro- 
genitor halo, due to the large amount of unbound particles. 
Besides the stringent requirement on the robustness of trac- 
ing, second order effects can lead to mass growth for satellite 
subhaloes. As we will show, both accretion and merger can 
continue to happen for satellites even within the virial radius 
of the host halo. Ignoring these process will under-estimate 
the mass of traced subhaloes. 

In this work we present a new tracing code (Hierar- 
chical Bound- Tracing, HBT hereafter) to identify subhaloes 
and construct the merger tree. The key to our code is a 
robust unbinding algorithm together with adaptive source- 
subhalo management. With these recipes we are able to walk 
through the cosmic age to capture every subhalo ever alive, 
as will be described in section (2] Our careful tracing of sub- 



haloes' hierarchy enables us to apply the unbinding algo- 
rithm recursively, naturally allowing satellite accretion and 
merger. The success of HBT is demonstrated in section [3] 
to have high completeness. Through comparison to a con- 
figuration space subhalo finder we also show that HBT sub- 
haloes are in general more robust and more extended. The 
results are summarized in section |4l HBT's implication for 
galaxy formation models and prospects are discussed in sec- 
tion [S] Two concordence LCDM simulations, both with cos- 
mological parameters Q.m = 0.3, SIa = 0.7, and as = 0.9, 
are used for the tests and comparisons in this work; one 
is a zoomed-in re-simulation of a IQ^^ h~^yiQ galax;y clus- 
ter with a particle mass = 1 x lO*/i~^M0 carried out 
with GADGET II (|Springef(|2005l )). and the other is a cos- 
mologica l simulation with a boxsize 100/i~^Mpc using 512 
particles (|jing fc Sutd[20o3 ) . 



Our algorithm also makes use of characteristics (ii) and 



(iii) and it turns out that the first characteristic is automat- 
ically satisfied. The difference between our algorithm and 
that of GIO is mostly technical, mainly in the time direction 
of tracing a subhalo: HBT proceeds from the earliest epoch, 
naturally finding the full hierarchy of subhaloes together 
with extracting the merger tree in only one walk through the 
cosmic time; while GIO tries to figure out the subhalo hier- 
archy level by level by subsequently revisiting earlier snap- 
shots. The MHT code dep ends on the 'MLAPM halo finder' 
(MHF; iGiU. Knebe. fc Gibson. 2004 ) to first identify haloes 
as well as subhaloes from the halo formation time. Then the 
haloes and subhaloes are tracked in subsequent outputs. In 
HBT the growth, merger, and stripping of haloes and sub- 
haloes are handled in a unified way starting from the earliest 
resolved haloes, and the full hierarchy of the subhalo merger 
tree is resolved. HBT also stands out in the sense that it is 
the first time that various systematic issues in a tracing al- 
gorithm have been investigated and the tracing results have 
been carefully assessed. 



2 ALGORITHM 

2.1 Overall tracing algorithm 

HBT starts the tracing of subhaloes from input halo cata- 
logues for a sequence of simulation outputs. Q In the cur- 
rent implementation we adopt the simple and widely-used 
Friends-of-Friends (FoF) halo finder to construct the in- 
put catalogues, because it has been demonstrated that viri- 
alised haloes and subhaloes are included in these FoF haloes 
(Springel et al. 2001), and also because it does not cut sub- 
haloes near halo bounda ries as, for example, an spherical- 
overdensity halo finder (|Lacev fc Cold (| 19941 )) would do. 
For our two LCDM simulations, the nominal linking-length 
6 = 0.2 is adopted. With halo catalogues in hand, the HBT 
algorithm can be summarized in one sentence: HBT builds 



^ Usually 60 snapshots for a LCDM simulation starting from 
an initial redshift where the first several haloes can be found, 
e.g. Zini ~ 20, is sufficient to achieve high completeness in the 
present-day subhalo catalogue. Interested readers can refer to Ap- 
pcndix [C3l whcre we investigate how the algorithm depends on the 
time resolution of simulation outputs. 
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and traverses the halo merger tree and finds the self-bound 
structure for every halo at every snapshot after its birth. 

Specifically, the subhalo identification breaks into two 
steps: first selecting candidate particles which contain the 
members of the self-bound structure, and then removing the 
unbound ones. We call the initial collection of the candidate 
particles, out of which the self-bound part gives rise to a 
subhalo, a source-subhalo (or source for short). In HBT we 
work with two types of subhaloes in each halo: a central 
subhalo which is the dominating subhalo within a host FoF 
halo, and satellite subhaloes which are the remaining ones 
if any. Starting from the highest redshift, haloes without 
progenitors are fed to an unbinding procedure as sources 
of their central subhaloes. At the next snapshot, the parti- 
cles of these source subhaloes are tracked to identify their 
host haloes. If more than one source is found to reside in 
the same host, a halo merger is identified The progeni- 
tor sources are then unbound to select the most massive 
self-bound structure as the central subhalo, and others as 
satellites. The sources for the central subhaloes are updated 
to be the current host halo excluding particles from any self- 
bound satellites, and unbound again to allow for growth of 
centrals. In this kind of tracing process, haloes without pro- 
genitors create new branches in the merger tree and give 
birth to new centrals, while mergers connect branches and 
transform centrals to satellites. 

Since the particles of satellite subhaloes are traced from 
their progenitors, this is equivalent to assuming that no ac- 
cretion of host halo's particles can happen for satellite sub- 
haloes. In Appendix |B] we explicitly test this assumption 
and show that local accretion has little effect on the tracing 
in the long term. 

We record the progenitor-descendant information as we 
proceed. This way we get the subhalo catalogues and the 
merger history of these subhaloes at the same time. 



2.2 Tracing Satellites Robustly and Efficiently 

As we have mentioned before, it is the most challenging part 
for a tracing algorithm to trace satellite subhaloes robustly 
and efficiently. This is a problem mostly because unbind- 
ing can be fastidious about the quality of source subhaloes. 
Either too big or too small a source subhalo can lead to fail- 
ure in unbinding. This can be understood as follows. The 
definition of boundness depends on the reference frame in 
which to calculate the kinetic energy and on the assembly 
of particles from which the potential energy is obtained. For 
a self-bound subhalo, the reference frame is defined to be at 
rest with respect to a certain "centre" of the subhalo (e.g, 
centre of mass) which removes the bulk motion of the sys- 
tem, and the potential energy is summed over all the par- 
ticles in the subhalo. Starting from a source subhalo with 
bound and unbound particles, one still seeks a centre which 
is roughly at rest with the final one, and uses this frame 
to calculate the binding energy. An overly large source sub- 
halo with too many unbound particles would probably give 



^ It also happens that a source's particles are distributed into 
multiple host haloes. In this case the source is split as described 
in Appendix ICll 



a centre which deviates too much from the underlying sub- 
halo, while an overly small one consisting of only a small 
portion of the real subhalo would give too shallow poten- 
tial, both yielding severely biased result or e ven completely 
missin g the subhalo. For example, as shown in lHavashi et al.l 
(|2003D an NFW halo truncated to a radius of < 0.77 times 
the scale radius is completely unbound. Thus a source sub- 
halo is expected to be a slightly enlarged assembly of parti- 
cles containing the final subhalo. In addition, the unbinding 
procedure is required to be robust, especially if the source 
is strongly contaminated. We make efforts in both aspects 
to improve the robustness of HBT. 

In HBT we implement a core-averaged unbinding algo- 
rithm designed to tolerate contamination. Starting from a 
source subhalo, unbound particles are removed iteratively 
till the bound mass converges. For each iteration, the ref- 
erence frame is chosen to be the centre of mass and bulk 
velocity of an inner-most core consisting of a certain frac- 
tion CoreFrac of the remaining particles with the lowest- 
potential. This provides a robust estimate of the underlying 
true reference frame in presence of contaminations. More 
detailed description of this algorithm can be found in Ap- 
pendix [X] 

Besides improving unbinding, our HBT code updates 
satellite sources adaptively, to keep the amount of contami- 
nating particles under control, while still maintaining a large 
enough reservoir of candidate particles. The idea is to reduce 
the size of the source conservatively every time a subhalo 
has been stripped to a fraction CoreFracO of the current 
source mass Msrc- To ensure that the new source with mass 
Msrc2 is still larger than the current subhalo, we replace the 
source with the subhalo's progenitor found at the time it 
was first stripped to a fraction \/CoreFracO of the original 
source mass. With the new source constructed this way, it 
is ready to see that the current subhalo automatically be- 
comes the new source in the next loop of source updates, 
saving the trouble to search for new sources upon every up- 
date. The source for a subhalo before infall into a host halo 
is simply its host halo with all the satellite particles masked 
out. After infall this kind of source updates start to take 
effect all the way along the subhalo's history. Note that the 
removal of unbound particles not only reduces contamina- 
tion, but also reduces the amount of calculation for unbind- 
ing. With source subhaloes updated this way, the unbinding 
procedure is always protected to work under the condition 
CoreFracO -Msrc < M^uh < V CoreFracO -Msrc, substantially 
reducing the amount of work for unbinding and making the 
tracing more robust. 

The CoreFrac parameter adopted by the unbinding pro- 
cedure is initially CoreFracO, and updated with CoreFrac = 
Msrc2/ Marc ' V CoreFracO every time we update the source. 
This is to ensure that the smallest subhalo out of Msrc 
does not use an over-large core when Msrc2 « Msrc ■ 
VCoreFracO due to discreteness in simulation output time. 
In the current implementation, we adopt CoreFracO = 0.25. 

In the development stage we also tried several other 
reference frames for unbinding, including: 

(i) centre of Mass (CoM) frame: Take the centre of mass 
of source particles as centre and the average velocity as the 
bulk motion. 

(ii) centre of Potential (CoP) frame: Take the potential 
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Figure 1. Robustness of various unbinding algorithms. Left: Subhalo mass evolution calculated wrhen the progenitor halo is used as a 
source vfithout any updates during tracing. The solid lines are results of eore-averaged unbinding, with different CoreFrac as given in 
the legend. The other three lines denote the results when three eommonly-used reference frames are adopted; the dotted line is for the 
centre of Mass(CoM) frame, the dot-dashed line for the MinPot frame, and the dashed line for the CoP frame. These reference frames 
are defined in the text. The thick grey line shows the halo-centric distance of the subhalo, in an arbitrary unit. Middle: The same as in 
the left column, but the subhalo at last tracing step is used as a source. Right: Subhalo mass evolution calculated when a self-adaptive 
source is used. The solid lines are the results of self-adaptive core-averaged unbinding with different CoreFracO parameters. Over-plotted 
are the results of the core-averaged unbinding with no source update (the dotted line, overlapping with the solid blue line in top panel; 
as in the left column) and with immediate updates (the dashed line; as in the middle column), with CoreFrac = 0.25. Top panels are for 
halo S51G86, while bottom panels for halo S43G11. 



weighted centre of mass as centre and the potential weighted 
average velocity as the bulk motion. 

(iii) Minimum Potential (MinPot) frame: Take the posi- 
tion of the particle with minimum potential energy as centre 
and average velocity of the source subhalo as the bulk mo- 
tion. We avoid using the velocity of the minimum potential 
particle considering that it may have large dispersion with 
respect to the bulk motion. 

For the source construction, we also tried two simple strate- 
gies. One is to always use the source from the infall with- 
out any updating (no- update), and the other is to update 



it aggressively by using the progenitor subhalo at snapshot 
n — 1 as the source for a subhalo at snapshot n (immediate- 
update). 

In Figure [T] we compare the performance of various 
combinations of these unbinding and source construction 
recipes, by applying them to two subhaloes in the clus- 
ter simulation. Top panels show the case for tracing halo 
S51G86 (named after snapshot ID and halo ID just before 
the infall, with the merger mass ratio 0.0016 and the initial 
orbital circularity 0.16). The left column is the "no- update" 
regime, where the source is taken to be the progenitor halo 
at the infall, while the middle column does "immediate up- 
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date" and uses the subhalo at the previous tracing step as 
the source. It can be seen that the CoM frame method and 
the MinPot frame method have similar performance; both 
methods fail to identify more than a half of the bound pop- 
ulation. The performance is worse at a smaller halo-centric 
distance where the tidal force is stronger or when the bound 
part is smaller hence there are more contaminating parti- 
cles. In general the core-averaged unbinding algorithm has 
the best performance, due to its ability to seek out a tight 
core to represent the majority of the bound particles. Still 
an over-estimation of the core size (e.g. CoreFrac — 0.5 after 
z — 0.2 with no update) can cause an obvious drop in subse- 
quent subhalo mass. As shown in the top left panel, the sub- 
halo gradually loses mass as it spirals into the centre of its 
host, reaches a minimum at pericentre passage near z — 0.6, 
and then re-gains mass as it moves out. Those regained par- 
ticles are also from the initial infalling halo, and we found 
that they are also mainly bound particles before the peri- 
centre passage. This rebinding of once unbound particles are 
suppressed in the top middle panel, where monotonic de- 
crease in subhalo mass is forced by keeping only bound par- 
ticles from last snapshot as the source. In an extreme case 
as shown in the bottom panels for halo S43Gll(with merger 
mass ratio 0.02 and initial orbital circularity 0.10), the re- 
binding process can increase the subhalo's mass by a factor 
of 2.3 after pericentre, and suppression of the re-capture can 
yield an under-estimation of subhalo mass by a factor of 2.8, 
comparing the left and middle panels. Right column shows 
the result when the adaptive update of source subhalo is 
applied. With the source size under control, the problem 
of an overly large core is avoided. Because Marc marks the 
maximum subhalo mass allowed, Adsrc < Msub/CoreFracO 
shows that the maximum factor by which a satellite subhalo 
is allowed to grow through rebinding is 1/CoreFracO. Thus 
a big CoreFracO would still suppress mass growth, as in the 
case of the bottom right panel with CoreFracO = 0.5. 

Figure [2] further shows the reason for the performances 
of different unbinding algorithms. At the time of infall, us- 
ing the reference frame of the final self-bound subhalo in 
the CoreFracO — 0.25 unbinding result as the standard one, 
we check the difference of the initial estimation of the ref- 
erence frames in each algorithm from that standard frame, 
and plot them normalized by the host halo virial radius Rvir 
and virial velocity Vvir ~ \/ GM^ir/Rvir- The CoM, Min- 
Pot, and CoP frames all have large velocity residuals with 
respect to the bound structure, although the MinPot frame 
has a negligible positional displacement However, the ini- 
tial CoM frame also has a large positional displacement. The 
assignment of these displaced frames would inevitably result 
in the artificial loss of bound mass in these three algorithms. 

The most-bound particle of an HBT subhalo is almost 
always located at the centre of mass position of the core, 
with a displacement within one softening length, reflecting 
the definition of our boundness. 
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Figure 2. The reference frames in different unbinding algorithms 
for halo S51G86. For each algorithm, we plot the difference of 
the initial estimation of the origin and the relative velocity of 
the reference frame from the final reference frame found for the 
subhalo in the HBT code, normalized by the halo virial radius 
and virial velocity. 



2.3 Allowing accretion and merger within satellite 
systems 

Although the dominating physical process in s atellite mass 
evolut ion is tidal mass loss, it has been shown in lSimha et al.l 
Hooi) that mass ac cretion and merger is not terminated 
for satellite galaxies. lAngulo et al.l |20o3) also find that the 
merger rate between satellite subhaloes and that between a 
satellite and a central subhalo are comparable for satellites 
smaller than 0.01 times the host halo mass, and that most 
of the satellite-satellite mergers happen between subhaloes 
which were once in central-satellite relation. These satel- 
lites and "satellite of satellite" s as well "satellite of satellite 
of satellite" s and so on define hierarchical satellite systems 
which coul d persist for a long period u p to many Gyrs as 
observed in lWhite. Cohn. &: SmitI (|2010l l. and these systems 
are places where satellite accretion and mergers can still hap- 
pen. 

In HBT we keep record of this "sat-of-sat" hierarchy. 
Each central subhalo is aware of its satellites, and each 
satellite is aware of its sat-of-sats which it had accumu- 
lated before its infall and which are still alive. There are 
occasions when a subhalo is ejected from a host , which is 
quite common a s shown in 'Lin. Jing. fc Li nI (|2003|) fsee also 

IOOTI:' ' 



Note that we also adopted the average velocity as reference for 
MinPot frame. 



Gill Knebe. fc Gibson 2005: Sales ct al, 20071: iLudlow et all 
20091) . In this case the ejected subhalo is also removed from 
the sat-of-sat list that contains it. Here we use the term 
"sat-o f-sat" instead of "sub-in-sub" (see.e.g. lSpringel et al.l 
(|2008t )') to emphasize that this relation between subhaloes is 
a historical or dynamical relation, which may not correspond 
to the spatial nesting relation in the latter case because of 
separation of orbits due to host halo's tidal force or mul- 
tiple satellite interactions. We allow the accretion of mass 
within each subhalo's sat-of-sats. This is achieved by imple- 
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meriting the unbinding procedure recursively: for a hierarchy 
of source subhaloes, starting from the highest level source 
subhalo, we feed the unbinding procedure with both par- 
ticles from the current source subhalo and particles which 
are removed from recursive unbinding of lower level source 
subhaloes|3 This means particles removed from a satellite 
can have a chance to be accreted into a higher level sub- 
halo. When a satellite subhalo dies with the majority of 
its particles accreted by a higher level satellite subhalo, a 
satellite-satellite merger happens. 

Ignoring the effect of satellite merger could result in loss 
of subhaloes during tracing. For example, a subhalo with 
7000 particles (or about 1/1000 the host halo mass) is miss- 
ing in the HBT result for the resimulation when satellite 
merger is switched off. This is triggered by a satellite major 
merger event between two satellites which were once in the 
central-satellite relation. After merger the particle kinetic 
energies are increased. If one still uses particles from only 
one source subhalo to do the unbinding, the potential would 
not be deep enough to bind the particles, resulting in the 
loss of the new subhalo. 

In Figure [3] we show the ratio of the total mass con- 
tained in each subhalo mass bin for the resimulated clus- 
ter, found by HBT with and without the satellite merger. 
This is equivalent to the ratio of the mass weighted subhalo 
mass function MsubdN/dln Maub{see e.g.. lGao et al.|[2004bl . 
for more details of the mass function). Although allowing 
merger among satellite subhaloes has almost no effect on 
small subhaloes, it can enhance the mass function at the 
high mass end by 20 percent, a fraction close to the frac- 
tional mass contribution from satellites to host haloes. The 
strong fluctuations at the high mass end reflect the occur- 
rence of significant satellite merger events. 



3 RESULTS AND COMPARISON 

To show that HBT has a superior tracing ability to recover 
subhaloes in high density regions and to compare with other 
subhalo finders, we have applied HBT t o one test case pro - 
vided by the "haloes Gone Mad" proiect (|Knebe et al.ll2011^ ■ 
to examine a subhalo's mass evolution. To show that HBT 
completely recovers the subhalo population, we compare the 
subhalo mass functions from HBT for our simulations with 
those from SUBFIND, as well as with a fitting formula. 
The comparison shows that HBT subhaloes are more mas- 
sive than those from SUBFIND by 10 to 20 percent, a fact 
that is also observed when comparing the size of subhaloes 
from the two codes. The reason for this difference is further 
revealed in the density profile, where SUBFIND subhaloes 
show sharp truncation near tidal radii while HBT subhaloes 
are more extended. We also show the one-to-one matching 
resuh between HBT and SUBFIND catalogues. HBT has 
also been compared with many ot her subhalo finders in the 
subhalo-finder comparison proiect (|Onions et al.ll20l3 ). and 
is found to have good performance. 



* To avoid adding satellite particles multiple times, we define the 
source of a central subhalo to be the host halo excluding all the 
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iources when recursive unbinding is 



Figure 3. The Effect of satellite merger. We plot the ratio of the 
total subhalo mass in each logarithmic mass bin between HBT 
results with {dM) and without (dMifSM) satellite merger for 
the resimulated cluster. The errors are propagated from Poisson 
errors in the total number of particles within each bin. Different 
colour lines with symbols are for different redshifts, while the 
horizontal dotted line is the 1 : 1 reference. 



3.1 Subhalo Mass Stripping History 

The test simulation we use from the "haloes Gone Mad" 
project is a head-on collision simulation of two NFW haloes 
with initial virial masses of IQ^^Mq and 1O^'*M0, corre- 
sponding to 10* and 10^ particles. The small halo is thrown 
right throug h the centre of the other halo. As shown in 
iKnebe et al.l (|201ll ). all halo finders based on only config- 
uration space information fail to find the subhalo at host 
centre, and even some finders based on position and veloc- 
ity phcise space information give an un-physical result near 
the centre. For clarity, we choose two representative config- 
uration based finders (AHF and SUBFIND) and two phase- 
space based finders (HSF and ROCKSTAR), and compare 
our HBT result in Figure|3]with those found by the four find- 
ers. It can be seen that the HBT subhalo is complete from 
the start, robust near the halo centre, and clean as it moves 
out of the central region. As just stated, the two configura- 
tion based finders fail to find the subhalo at the very central 
region of the host halo. HSF gives an unphysically high sub- 
halo mass in the beginning (higher than the halo mass before 
the infall), which may be attributed to the inclusion of some 
local particles as discussed in section [B] AHF gives a much 
higher mass when the subhalo just passes the centre, which 
is not physical because, as can be seen from the right panel, 
the maximum circular velocity Vmax well exceeds that of the 
progenitor halo. A significant fluctuation in the Vmax history 
is also observed for the subhalo found by ROCKSTAR near 
the host centre, indicating some ambiguity in capturing the 
bound mass by the finder. Furthermore, there is no obvious 
reason why the subhalo's mass should peak at 0.4 times the 
host virial radius as given by ROCKSTAR. In general HBT 



Hierarchical Bound-Tracing Algorithm 7 



HBT 
AHF 

SUBFIND - 
HSF 

ROCKSTAR 




200- 



150- 




! 100- 



50- 



HBT 
-AHF 

-SUBFIND 
-HSF 

-ROCKSTAR 



1.4| , , 1 250r 

1.2 
1 

.0.8 
0.6 
0.4 
0.2 


-2-1 1 2 

D/R 

virMost 

Figure 4. Mass stripping liistory of a subhalo falling through the centre of a host halo as found by HBT and the other four representative 
subhalo finders. Left: Bound mass of the subhalo normalized by its virial mass before infall, as a function of its distance from host centre, 
normalized by the virial radius of the host halo. The subhalo moves toward the positive direction of the axis. Right: Evolution of the 
maximum circular velocity of the subhalo. The circles on the ROCKSTAR line mark the timesteps of the simulation outputs. 



has comparable performance with HSF, while still has the 
advantage that all the HBT particles are strictly physically 
associated with the subhalo by design. 

3.2 Subhalo mass function 

In Figure [5] we show the subhalo mass function found by 
HBT and SUBFIND for the resimulated cluster at z=0. 
For reference we also plot the power - law fitting formula 
jGao et al.ll2004bl: ISpringel et aLlliooi : lAngulo etaLll2009l : 
iGiocoh et aLlbOld l. 



dN 



where Msub is subhalo mass and Mhost is the host halo virial 
mass. The normalization A'o depends on the definition of the 
virial radius within which subhaloes are counted. For the 
virial relation predicted by spherical collapse model(see e.g. 
iBrvan fc Normajilll998l . for a fitting formula), GIO found the 
normalization to be A^o = lO"^ ''^(/i^^M0)"° ^ We adopt 
the same virial definition throughout this paper, with the 
average density within the virial radius defined to be 101 
times the critical density at = 0, and will compare our 
subhalo mass function with the GIO result. It can be seen 
that the mass functions are consistent with the fitting for- 
mula from GIO, with HBT's subhaloes being more massive 
than the result of SUBFIND by 10 to 20 percent. 

3.3 Size of Subhaloes 

The size of a subhalo is trimmed by the tidal force from 
the host halo, and can be estimated by a tidal radius which 
is defined as the radius of the satellite subhalo at which 
its self-gravity equals the tidal force of the host halo. In the 
limit Rtidai « D where D is the halo-centric distance of the 
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Figure 5. Subhalo mass function at z=0 for the cluster resim- 
ulation. Top panel shows mass weighted subhalo mass function 
as found by HBT and SUBFIND. The black solid line is the fit- 
ting formula of Eq uation 13.11 with the normalization given by 
iGiocoli et al.l l|201Cll V Bottom panel plots the ratio of the total 
subhalo mass contained in each mass bin from the two codes, or 
equivalcntly the ratio of the mass-weighted subhalo mass func- 
tions. The horizontal dashed line is the 1:1 reference. 



subhalo, the tidal radius is estimated as (|Binnev fc Tremaind 
1 19871 : iTormen. Diaferio. fc Sveill 19981 ): 



(2- 
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Figure 6. Size of subhaloes in the 100 most massive lialocs in 
the cosmological simulation as found by HBT and SUBFIND. 
Reg and Req,Q.02 are the radii where the subhalo density equals 
to 1 and 0.02 times the background density. Rudal is the subhalo 
tidal radius. We use only FoF haloes containing more than 1000 
particles and subhaloes containing more than 100 particles. The 
red solid lino in each panel marks the median value of Req (or 
^eij,0.02 according to y-axis). In the lower left panel we overplot 
the median line of -Reg, 0.02 of SUBFIND result in green dashed 
line. 



For singular isothermal density profile, it is easy to 
check from Equation (|3.2p that the local density of the host 
halo equals to that of the satellite at tidal radius. Thus 
a radius of equality R^q can be defined at which the den- 
sity of the satellite is the same as the local density of the 
host. For realistic density profiles, the tidal r adius may dif- 
fer fro m R^q by a factor of order unity. In ISpringel et al.l 
(|2008l ) they found that the subhalo tidal radius equals to a 
radius Jie9,o.02 where the subhalo's density falls below 0.02 
times the local host density. We compare the tidal radius 
to Req and iie9,o.02 for both HBT and SUBFIND results in 
Figure El using the 100 most massive FoF haloes in the cos- 
mological simulation. To avoid resolution effect, we use only 
FoF groups with more than 1000 particles and subhaloes 
with more than 100 particles. It can be seen that although 
the Re.qfi.02 is very close to Rudai for SUBFIND result, it 
is in general 20 ~ 30 percent larger in the HBT case. But 
for Req the two codes give consistent results. This reflects a 
puff-up in the outer region of subhaloes in HBT result, or 
rather a sharp cut-off in the SUBFIND result, due to particle 
division using density saddle points. We show this explicitly 
in the next subsection. 



3.4 Density Profile 

There is almost no difference in the density profile between 
central subhaloes found by HBT and SUBFIND. But for 
satellite subhaloes, especially when they reside in the central 
region of host haloes, HBT gives a much more extended pro- 
file. Figure [7] compares the density profiles of the most mas- 
sive satellite of the cluster found with the two codes. This 
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Figure 7. Density profiles for the most massive satellite sub- 
halo of the resimulated cluster as found by HBT and SUBFIND. 
Red lines are for particles which are bound to the subhalo, green 
lines are the total density of all the other particles. The vertical 
black lines mark the tidal radius of the subhalo. Solid lines are 
the HBT result while dashed ones are the SUBFIND result. The 
vertical dotted line on the left marks the smoothing length of the 
simulation. 



subhalo is resolved to have lO" particles or 1/10 the virial 
mass of the cluster. It is located at 0.37?«ir from the cluster 
centre. As expected, a sharp cut is observed for the SUB- 
FIND profile near its tidal radius, while HBT gives a subhalo 
mass which is 6.5 times as large. This is due to HBT's prior 
knowledge of source particles beyond the tidal radius while 
SUBFIND can only identify those particles within tidal ra- 
dius, as already seen in section [331 



3.5 Cross Match between HBT and SUBFIND 

To see if the two codes find the same set of subhaloes, we try 
to match the members in the two catalogues. Given one tar- 
get subhalo in hand, we search host subhaloes for its member 
particles in the other catalogue. For each host subhalo, we 
calculate a boundness weighted summation of the number of 
matched particles, with the most-bound particle in the tar- 
get subhalo having the largest weight. The host subhalo with 
the biggest summation is selected as the target's correspon- 
dence in the other catalogue. We do the match from HBT 
to SUBFIND and vice versa. We classify the subhaloes into 
three categories according to the two matches: bilaterally 
matched subhalo when the subhalo and its correspondence 
are matched to each other; unilateral subhalo when the cor- 
respondence cannot be matched back to the subhalo; and 
un-matched subhalo when no correspondence can be found 
in the other catalogue. 

The match result for the galaxy cluster is shown in Fig- 
ure [HI The majority of subhaloes are bilaterally matched. 
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accounting for 82% HBT subhaloes and 97% SUBFIND sub- 
haloes. 90% of the bilateral SUBFIND subhaloes have more 
than 95% of their particles shared with their HBT corre- 
spondences. The median line points out 10 to 20 percent 
increase in mass for HBT subhaloes compared to their SUB- 
FIND correspondences. Only less than 1 percent subhaloes 
belong to the un-matched category, with their masses being 
lower than 100 particles. 

Unilateral subhaloes are discussed extensively in Ap- 
pendix [D] Here we only focus on one particularly interest- 
ing case. When two haloes with comparable masses merge, 
it is likely that the resulting halo contains two subhaloes 
with comparable peak densities, making it difficult to dis- 
tinguish which is superior and which is subordinate without 
appealing to progenitor information. Especially when these 
two subhaloes are close to each other forming a high density 
environment, the difficulty of resolving and weighing them 
are elevated. We give an example of such a binary system 
in Figure [9] from the cosmological simulation. Two satellite 
subhaloes sit at the centre of the image resulting from a 1:2 
merger of two progenitor haloes. This binary satellite system 
is found near the boundary of the host halo and thus closely 
resembles an individual halo. Most of the surrounding parti- 
cles are bound to both of the density peaks. SUBFIND takes 
the peak from the smaller halo which has higher central den- 
sity and associates the surrounding particles with it to make 
a parent subhalo, leaving the core from the bigger halo as 
its sub-in-sub, while HBT gives a reversed sub-in-sub hierar- 
chy which is consistent with the merger hierarchy. When we 
match these two subhaloes from SUBFIND to HBT, obvi- 
ously the small subhalo in SUBFIND is matched to the big 
one in HBT. And because the surrounding mass outweighs 
the core in the big subhalo of SUBFIND, this one is also 
matched to the big subhalo in HBT. For the same reason, 
the two HBT subhaloes are also matched to the big SUB- 
FIND subhalo. This may be of particular importance when 
constructing merger histories. A false switch between two 
branches of a tree would happen as a result of the switch 
between the sub-in-sub relation of the binary. With the help 
of subhaloes' merger history, HBT is able to give more phys- 
ical weights to the dominance of the two cores and avoid this 
kind of faulty links in the merger tree. In Figure [SI these bi- 
nary subhaloes would produce two adjacent perpendicular 
lines connecting near the 1:1 line if the surrounding parti- 
cles could outweigh one core, or a single line if one subhalo is 
missing. A dedicated study of how major-mergers challenge 
different subhalo finders can be found in Behroozi et al.(in 
prep). 



4 SUMMARY 

In this work, we have managed to develop a tracing algo- 
rithm to produce a complete and physically-motivated sub- 
halo catalogue. In our HBT algorithm, subhaloes are di- 
vided into two types where central subhaloes grow via ac- 
cretion from the host halo and from satellites inside the 
same host while satellite subhaloes can only accrete from 
within their satellite-of-satellites. We keep a full record of 
subhaloes' merger hierarchy and apply the unbinding algo- 
rithm hierarchically to enable satellite accretion and merger. 
While omitting satellite accretion can result in 20 percent 




^SUBFIND /particles 

Figure 8. Cross match between HBT and SUBFIND subhaloes 
for the cluster. Each subhalo is plotted according to its mass and 
the mass of its correspondence. Blue dots are bilaterally matched 
subhaloes (see text for definition). The solid cyan line gives the 
median masses of the bilateral subhaloes. Red circles are unilat- 
eral HBT subhaloes and green squares are unilateral SUBFIND 
subhaloes, except for those on the axes which are unmatched 
ones. For each unilateral match, we also draw a line connecting 
that subhalo to its correspondence. The point at the top-right 
corner is the central subhalo. 



loss in mass for massive satellites, individual tests of subhalo 
mass loss history and statistical comparison with SUBFIND 
which is based on local density search show that local accre- 
tion of background particles have only temporary effect and 
is negligible for satellite subhaloes. 

Because of the large dynamic range in the evolution of 
subhaloes' masses, a robust unbinding algorithm is required 
to trace subhaloes long enough. We have proposed core- 
averaged unbinding algorithm together with self-adaptive 
update of source subhaloes to accomplish this job. The core- 
averaged unbinding algorithm succeeds in its ability to give 
quick and accurate estimation of the centre and bulk veloc- 
ity of the subhalo, allowing safe removal of unbound par- 
ticles. Self-adaptive update of source subhaloes keeps the 
dynamic range in a safe region for unbinding while allowing 
rebinding of stripped particles. We find that successive trac- 
ing of satellites without allowing rebinding can cause up to 
a factor of 3 suppression in the subhalo mass at apocentre 
passage. Application of HBT to a contrived simulation also 
demonstrates its superb performance in robustly recovering 
the subhalo's stripping history. 

Since our algorithm utilizes historically constructed 
source subhaloes rather than locally collected sources using 
density and position thresholds, we can resolve satellite sub- 
haloes well even in the central high density region of haloes 
while SUBFIND may have severe spatial truncation for sub- 
haloes in halo centre. Even for small satellites our subhalo 
mass function is about 15 percent higher than SUBFIND. 
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Figure 9. An example of sub-in-sub switch. These are two sub- 
haloes forming a binary system near the boundary of the host 
halo. The lower-left panel shows the projected squared density of 
this system. The red dashed circle in the lower-left panel marks 
the virial radius of the system. The upper-right panel shows the 
contour map of the density field for the very inner region. Red 
symbols mark the most-bound particles of two HBT subhaloes; 
green symbols mark those of SUBFIND subhaloes. Pentagrams 
represent the superior subhaloes and crosses stand for subordi- 
nate subhaloes. The masses for the upper and lower subhaloes 
in the image are 94468(488) and 17541(111884) particles in the 
HBT(SUBFIND) case. 

The size of subhaloes found by HBT extends 20 to 30 per- 
cent larger than SUBFIND characterized by Reqo.o2- 

Because we do not need to do density interpolation 
or spatial searching to construct source subhaloes, HBT 
runs fast. The HBT code is written in C and has been 
fully OpenMP parallelized both for cosmological simulations 
where the parallelization is on halo level and for high res- 
olution resimulations where the parallelization is done on 
particle level; a hybrid MPI/OpenMP version is under de- 
velopment and will be available in the near future. It would 
also be straight forward and efficient to integrate HBT to a 
simulation code for on-the-fiy high resolution subhalo find- 
ing and merger tree output. In that case HBT would also 
not be limited by the number of snapshots in the simulation 
output. 



5 DISCUSSION 

Much of the treatment in HBT is more physical than math- 
ematical. Its complexity lies in the physical process involved 
and its success reflects its correct understanding of the sub- 
halo's evolution. This code is obviously not applicable to hot 
dark matter simulations, where structure formation is top- 
down rather than bottom-up. We would also expect HBT 
to have some difficulty when applied to warm dark mat- 
ter simulations, w here fragmentation is a vital p rocess in 
forming structures ()Bode. Ostriker. fc Turokl l2001f l. though 
our splitting algorithm would alleviate the problem. In fact 



HBT can be regarded as a low-level semi-analytic model for 
cold dark matter subhaloes. Its high resolution and phys- 
ical particle partitioning guarantees that the merger trees 
built by HBT have much fewer lost nodes or false links. 
This makes HBT an ideal choice for building merger trees 
for Semi- Analytic Models of galaxy formation. Besides, the 
tracking nature makes HBT easily extensible to find not only 
bound subhalos, but also more general structures such as 
streams. A combination o f HBT's tracking ability with the 
STructure Finder's (STF: lElahi. Thacker. fc Widrowir201ll ) 
stream identification algorithm is under development to 
trace streams (Elahi et.al. 2012, in prep). 

One example revealing HBT's advantage for galaxy for- 
mation models is the "sub-in-sub switch" case shown in Fig- 
ure [9] and discussed in section 13.51 When two subhaloes of 
comparable mass form close pairs, HBT is superior in its 
ability to partition the surrounding particles according to 
their origins. This is important because the surrounding par- 
ticles, while being co-bound by the binary-subhalo and can 
be assigned arbitrarily to either member, play a vital role 
in establishing links to subhaloes' progenitors and descen- 
dants, especially when they outweigh the cores of the bi- 
nary. HBT's way of partition would guarantee that these 
surrounding mass obey the progenitor-descendant relation, 
while a partition without knowing the origin of these parti- 
cles can easily lead to an incorrect link of the cores to the 
progenitor haloes. 

One future direction for the improvement of HBT code, 
and other subhalo finders as well, would be to find a bet- 
ter definition for subhaloes. Even with our ability to clearly 
construct source subhaloes near halo centre, subhaloes can 
still appear over-stripped near pericentre after which some 
stripped mass get rebound, posing a remaining "resolution" 
problem. We argue this remaining "resolution" problem to 
be due to the operational definition of subhaloes as instan- 
taneously self-bound structures, which is currently adopted 
by almost every subhalo finder. This definition only guar- 
antees coherence for a structure in isolation and in the col- 
lisionless limit of its member particles. Putting aside the 
structural evolution of the subhalo itself, there is still an 
evolving background which lowers the gravitational poten- 
tial of t he subhalo while a t the same time produces tidal 
force. In IShaw et al.l (|2007h both the tidal force from the 
host and the gravitational potential of stripped particles are 
considered to improve the coherency of subhalo definition. 
They find that the contribution from the tidal force and 
that from the lowered potential due to background particles 
roughly offset each other, with a slight increase for the mass 
of subhaloes in the inner region of a halo. 

It should be noted that tidal force is not always disrup- 
tive. This could easily be illustrated in the following simple 
picture. Consider two test particles with infinitesimal mass 
moving on the same elliptical orbit around a central mass 
M, starting from apocentre with a small time delay dt be- 
tween them. The relative velocity is then dlf = ^^^l^dt. 
Hence the internal kinetic energy of the system varies as 
dv^ oc r~*, reaching a minimum at apocentre and maxi- 
mum at pericentre. Because the internal gravity of the two 
test particles can be ignored, the change in the internal ki- 
netic energy is solely from work done by tidal force. It is 
ready to see that the tidal work helps to reduce the relative 
kinetic energy of the two particles from pericentre to apoc- 
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entre. As a result, particles previously marked as unbound 
to a subhalo can become bound again as the subhalo moves 
to the outer region of its host. This has been seen in our case 
study of the evolution of individual subhaloes in section [2]2j 
where strong rebinding is observed as satellites move out 

from the centres of haloes. 

We comment that in IShaw et al.1 l|200'if ) even though 
tidal energy has been included in the unbinding procedure, 
they remove velocity outliers before adding tidal energy, still 
failing to allow these high velocity particles to be re-captured 
through tidal work. 



ACKNOWLEDGEMENTS 

We thank Volker Springel for making the SUBFIND code 
available, and the anonymous referees for lots of help- 
ful comments. This work in Shanghai is supported by 
NSFC (10821302, 10878001, 11033006), by 973 Program 
(NO.2007CB815402), and by the CAS/SAFEA International 
Partnership Program for Creative Research Teams (KJCX2- 
YW-T23). HYW is Supported by NSFC 11073017. During 
the final stage of this work, JXH is supported by the Eu- 
ropean Commissions Framework Programme 7,through the 
Marie Curie Initial Training Network Cosmo-Comp (PITN- 
GA-2009-238356). 



REFERENCES 

Angulo R. E., Lacey C. G., Baugh C. M., Frenk C. S., 2009, 

MNRAS, 399, 983 
Barnes J., Hut P., 1986, Natur, 324, 446 
Barnes J. E., Hut P., 1989, ApJS, 70, 389 
Baugh C. M., 2006, RPPh, 69, 3101 

Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2011, 

arXiv:1110.4372 
Benson A. J., 2010, PhR, 495, 33 
Binney J., Tremaine S., 1987, gady.book. 
Bode P., Ostriker J. P., Turok N., 2001, ApJ, 556, 93 
Bryan G. L., Norman M. L., 1998, ApJ, 495, 80 
Davis M., Efstathiou C, Frenk C. S., White S. D. M., 1985, 

ApJ, 292, 371 

Diemand J., Moore B., Stadel J., 2004, MNRAS, 352, 535 
Efstathiou C, Davis M., White S. D. M., Frenk C. S., 1985, 
ApJS, 57, 241 

Elahi P. J., Thacker R. J., Widrow L. M., 2011, MNRAS, 
418, 320 

Gao L., De Lucia G., White S. D. M., Jenkins A., 2004, 

MNRAS, 352, LI 
Gao L., White S. D. M., Jenkins A., Stoehr F., Springel 

v., 2004, MNRAS, 355, 819 
Gao L., Frenk C. S., Boylan-Kolchin M., Jenkins A., 

Springel V., White S. D. M., 2011, MNRAS, 410, 2309 
Ghigna S., Moore B., Governato F., Lake G., Quinn T., 

Stadel J., 2000, ApJ, 544, 616 
Ghigna S., Moore B., Governato F., Lake G., Quinn T., 

Stadel J., 1998, MNRAS, 300, 146 
GiU S. P. D., Knebe A., Gibson B. K., 2004, MNRAS, 351, 

399 

GiU S. P. D., Knebe A., Gibson B. K., 2005, MNRAS, 356, 
1327 



Giocoli C, Tormen G., Sheth R. K., van den Bosch F. C, 

2010, MNRAS, 404, 502 
Giocoh C, Tormen G., van den Bosch F. C, 2008, MNRAS, 

386, 2135 

Hayashi E., Navarro J. F., Taylor J. E., Stadel J., Quinn 

T., 2003, ApJ, 584, 541 
Jing Y. P., 2002, MNRAS, 335, L89 
Jing Y. P., Suto Y., 2002, ApJ, 574, 538 
King I., 1962, AJ, 67, 471 

Klypin A., Gottlober S., Kravtsov A. V., Khokhlov A. M., 

1999, ApJ, 516, 530 
Knebe A., et al., 2011, MNRAS, 415, 2293 
KnoUmann S. R., Knebe A., 2009, ApJS, 182, 608 
Lacey C, Cole S., 1994, MNRAS, 271, 676 
Lin W. P., Jing Y. P., Lin L., 2003, MNRAS, 344, 1327 
Ludlow A. D., Navarro J. F., Springel V., Jenkins A., Frenk 

C. S., Helmi A., 2009, ApJ, 692, 931 
Maciejewski M., Colombi S., Springel V., Alard C, Bouchet 

F. R., 2009, MNRAS, 396, 1329 
Moore B., Governato F., Quinn T., Stadel J., Lake G., 

1998, ApJ, 499, L5 
Moore B., Ghigna S., Governato F., Lake G., Quinn T., 

Stadel J., Tozzi P., 1999, ApJ, 524, L19 
Muldrew S. I., Pearce F. R., Power C, 2011, MNRAS, 410, 

2617 

Nagai D., Kravtsov A. V., 2005, ApJ, 618, 557 

Onions J., et al., 2012, MNRAS, 423, 1200 

Sales L. V., Navarro J. F., Abadi M. G., Steinmetz M., 

2007, MNRAS, 379, 1475 
Shaw L. D., Weller J., Ostriker J. P., Bode P., 2007, ApJ, 

659, 1082 

Simha V., Weinberg D. H., Dave R., Gnedin O. Y., Katz 

N., Keres D., 2009, MNRAS, 399, 650 
Skibba R. A., van den Bosch F. C, Yang X., More S., Mo 

H., Fontanot F., 2011, MNRAS, 410, 417 
Springel V., et al., 2008, MNRAS, 391, 1685 
Springel V., Yoshida N., White S. D. M., 2001, NewA, 6, 

79 

Springel V., 2005, MNRAS, 364, 1105 

Springel V., White S. D. M., Tormen G., Kauffmann G., 

2001, MNRAS, 328, 726 
Tormen G., 1997, MNRAS, 290, 411 

Tormen G., Diaferio A., Syer D., 1998, MNRAS, 299, 728 
Tormen G., Moscardini L., Yoshida N., 2004, MNRAS, 350, 
1397 

White M., Cohn J. D., Smit R., 2010, MNRAS, 408, 1818 



12 J. Han et al. 



APPENDIX A: UNBINDING ALGORITHM 

The core-averaged unbinding in HBT is done in the foUowing 
steps. 

(i) Starting from an initial assemble of A*" particles of a 
source subhalo. 

(ii) We calculate the gravitational potential -0 for each 
particle from the contribution of all the oth er A'^ — 1 parti- 



cles. We us e the Barnes-Hut tree algorithm ( Barnes fc Hut 
19861 . Il989ll as im plemented in GADGET lfSpringel et al 



200ll : ISpringeill2005h to calculate the potential. We then find 
a fraction (set by a parameter CoreFrac) of particles with 
the lowest potential, and call it as the lowest potential core. 

(iii) We calculate the kinetic energy K of each particle, 
including Hubble flow, with respect to the average velocity 
and centre of mass of the lowest potential core. Then we 
remove any particles with a positive total energy K + tp > 0. 

(iv) If the remaining number of particles A^;, is below the 
user desired mass limit NBoundMin, the unbinding procedure 
stops with no subhalo found. Otherwise, if A'';, agrees with N 
within some accuracy MassPrecision, unbinding stops with 
these Nb particles as a subhalo. In the other case, take the 
remaining Nt particles as an initial assemble of particles and 
repeat the above steps 2-4. 

The parameter MassPrecision is introduced to im- 
prove the efficiency of unbinding, which is the most time- 
consuming part of the entire code. Strictly speaking the def- 
inition of self-boundness requires the iterative removal of 
un-bound particles to be 100 percent complete, i.e, the it- 
eration can only stop when the number of bound particles 
Nbound equals to the number of remaining particles A". How- 
ever, as we have discussed in section [5] self-boundness is not 
a perfect definition to identify particles associated with a 
subhalo. Practically, it would be enough to just give an es- 
timation of the bound structure which is fairly close to the 
final self-bound part. 

We investigated the speed of convergence for several 
randomly selected subhaloes and find that the bound mass 
converges quickly with the number of iterations towards the 
final self-bound mass, as shown in Figure [All After several 
iterations, the majority of unbound particles have already 
been removed. Further iterations only refine the bound mass 
to a higher accuracy, with smaller and smaller change in 
mass as the iteration goes on. An estimation of the precision 
of bound mass can be obtained as the ratio of bound mass at 
two subsequent iterations. By stopping the iteration at some 
precision threshold MassPrecision, a lot of iteration steps 
can be saved while keeping the mass estimation accurate to 
the pre-set accuracy. We adopt MassPrecision = 0.995 in 
the current implementation of HBT. 
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Figure Al. Typical convergence curve of a sublialo during un- 
binding. Top: The remaining bound mass after each iteration; 
Bottom: Estimated precision of the bound mass after each itera- 
tion. 



radius (defined as if the subhalo is a halo). Except in the 
first several snapshots (and the first apocentre passage in 
the S43G11 case) when the subhalo is around the boundary 
of its host, no accretion of local background mass is found. 



APPENDIX B: POSSIBLE LOCAL ACCRETION 

The physical assumption under tracing algorithms is that 
the mass accretion of satellite subhaloes can be ignored 
within a host halo, so that one needs only to follow the 
remnants of infalled haloes. In Figure IBll we give a direct 
test of this assumption. After extracting a self-bound sub- 
halo from its adaptive source, we search for and add any 
additional bound particles located within 2 times its virial 



APPENDIX C: MAXIMIZING THE 
RESOLUTION OF HALOES AND SUBHALOES 

CI Splitting algorithm 

There are several conditions under which one may need to 
split a source subhalo into parts. One case is that a source 
subhalo may contain additional substructures which are un- 
resolved given the time and spatial resolution of the simu- 
lation outputs, or from a casual link of two haloes as one. 
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Figure Bl. Effect of local background particles vs. infalled 
source subhalo. The solid line shows the bound mass from 
the infalled halo, as done in the right panel of Figure [T] with 
CoreFracO = 0.25. The dashed line shows the mass after adding 
bound local particles. The masses are normalized by the sub- 
halo's mass at infall. Thick grey and thin black lines are for haloes 
S43G11 and S51G86 respectively. 



These unresolved substructures may separate from the host- 
ing source subhalo later in orbits, due to orbital parameter 
difference or tidal forces of nearby structure. Also a strong 
tidal force can break up a nearby smooth halo. In these cases 
a source subhalo may have multiple descendants hosted by 
several haloes. We implement the following splitting algo- 
rithm to take this into account: 

For each source subhalo, we partition its particles into 
fragments according to their current host haloes. Fragments 
with particle numbers smaller than NSrcMin are discarded. 
Here NSrcMin is a parameter controlling the mass resolution 
of source subhaloes. We choose this to be equal to NBoundMin 
in our implementation. If all the fragments are discarded, 
discard this source subhalo. The subhalo out of the most 
massive fragment is taken as the "main descendant" of this 
source subhalo. Other remaining fragments are called splin- 
ter source subhaloes, the subhaloes directly out of which are 
named splinter subhaloes. We keep a record of both the main 
descendant and the splinter information, yielding a merger 
tree which can have "downward branches". Note that be- 
cause we do not appeal to spatial clustering within haloes 
to identify subhaloes, our tracing algorithm can only re- 
solve splitting when it happens spanning over several haloes. 
Within one halo even if a subhalo breaks into two obvious 
part, as long as they are still bound as a whole we have no 
way to split them in HBT. 

The mass function of the biggest splinters (excluding 
the main descendant) in each splitting event decays approx- 



imately as 



much steeper than the slope of 
- oc m""'®. It's ready to see that 



d In m 

subhalo mass function 
these splinters are primarily small subhaloes near the reso- 
lution limit. 



C2 Quasi-halo treatment 

We take those subhaloes for which we cannot find their host 
FoF haloes as being in a "background halo" when assigning 
hosts in tree constructing algorithm as well as when splitting 
source subhaloes. This way we can find subhaloes which do 
not have a known host halo. Their source subhaloes can be 
regarded as "quasi-haloes" which are not identified by the 
halo-finder, but they still have concentrated structure and 
progenitor-descendant link between neighboring snapshots. 
These quasi-haloes are haloes fluctuating around the resolu- 
tion limit. The biggest quasi-haloes identified have about 50 
particles when haloes are filtered by a mass limit of 10 parti- 
cles in our implementation. Although they could be split or 
ejected subhaloes, those who can grow to a meaningful size 
are almost all haloes at their earliest stage of growth, i.e, 
haloes which have just become resolvable. The introduction 
of quasi-haloes is intended to recover the lost nodes in the 
merger tree due to halo-finder pitfalls, but has little effect 
on the tree even if omitted because they come as primarily 
early fluctuations. 



C3 Time resolution dependence 

In our tracing algorithm, one problem of concern is whether 
the time resolution of simulation outputs would affect the 
tracing result. The effect of time resolution is two-fold: time 
resolution translates into resolution of halo growth history; it 
also changes the dynamic range for unbinding thus affecting 
its robustness. 

First, there could be un-resolved satellites as well as 
un-resolved growth of satellites. Obviously those haloes born 
between two subsequent outputs and immediately becoming 
satellite subhaloes in the second snapshot are not recorded 
in halo catalogues and their descendants as independent sub- 
haloes are missed. However, since these haloes get stripped 
right after they become discernible, they do not have time 
to grow and are always fairly small haloes. Also the bet- 
ter the time resolution and the smaller the halo mass limit, 
the smaller are these unresolved haloes. For those haloes 
which merge between two subsequent snapshots, the growth 
of satellite haloes after the first snapshot and before merger 
are not captured either. Thus one would expect a decrease 
in subhalo population with decreasing time resolution. 

To quantify this effect, we compare the subhalo mass 
function M^^^_^dN / dluMsub identified with different time 
resolution. Our cosmological simulation has 60 outputs from 
z — 15.02 to z = equally spaced in log-space of the scale 
factor. We dilute these outputs by keeping one snapshot af- 
ter skipping every a few number of snapshots. Subhaloes 
are identified using these diluted snapshots with HBT. It is 
found that the diluted subhalo mass function has the same 
shape as the best resolved one, with the ratio of them well fit- 
ted by a constant and depending weakly on host halo mass, 
as long as the time step used for tracing is not too large. 
In Figure ICll we examine the relative amplitude of the di- 
luted subhalo mass function to the best resolved one (with 
Alna = 0.0468) at five different redshift for the cosmolog- 
ical simulation. The relative amplitude A as a function of 
time resolution can be well fitted by a Gaussian function 
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Figure CI. The Effect of time resolution on subhalo abundance. 
Data points are the relative amplitudes of the diluted subhalo 
mass functions to the best resolved one, as a function of time 
resolution. We show the result at five different redshifts, together 
with a Gaussian function fitting at each redshift. Error bars in 
the figure are propagated from Poisson errors on the subhalo mass 
functions. 



with parameters Ao and b, where A In a is the difference in 
In a between subsequent snapshots with a being the scale 
factor. The parameter Aq represents the relative amphtudes 
extrapolated to infinitesimal time resolution. The fitted Aq 
are all fairly close to unity, increasing from 1.001 ± 0.001 
to 1.003 ± 0.001 from 2 = to 2 = 2.1, indicating that 
the best resolved mass function is over 99.7 ~ 99.9 percent 
complete. The parameter b depends on the scale factor as 
b = 1.29a + 0.51. Equation ICll provides a handy reference 
in assessing the completeness A/Aq for subhalo populations 
from HBT, although the exact completeness for a particular 
application would also depend on the mass resolution of the 
halo and subhalo catalogues. 

The other aspect of the time resolution problem is 
whether the unbinding routine would be dependent on the 
time resolution adopted. We still use the halo S51G86 stud- 
ied before to test this. As shown in Figure [C2] applying the 
adaptive unbinding algorithm with different time resolution 
has no noticeable effect on the resulting self-bound mass. 
This again demonstrates the robustness of our core-averaged 
unbinding algorithm in conjugation with our adaptive source 
management. 



APPENDIX D: 
SUBHALOES 



UNILATERALLY-MATCHED 



Unilateral subhaloes in Figure [8] come in several situations. 
First there are subhaloes found by one finder but missed 
by the other, and thus they appear as part of a larger sub- 
halo, either central or satellite. These are refiected as vertical 
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Figure C2. Effect of time resolution on the self-bound mass. We 
apply the code skipping different number of snapshots at every 
tracing step, to test the stability of our unbinding algorithm with 
respect to time resolution. In the top panel, the solid line denotes 
for the bound mass when using all the snapshots available, while 
different symbols give the results of skipping different numbers of 
snapshots, labelled as different steps in the scale factor a. In the 
lower panel, we show the relative difference with respect to the 
highest resolution case. 



lines connecting up or horizontal lines connecting right to bi- 
lateral subhaloes. Most unilateral subhaloes belong to this 
case, with the majority having their correspondences being 
the central. For those unilateral HBT subhaloes linked to 
a satellite, most of them are in fact overlapping with their 
correspondence, having a separation between their centres 
being the order of the smoothing length of dark matter parti- 
cles. In fact if they can keep this overlapping state, these two 
subhaloes should be treated as one, as merger between them 
has completed. However, because overlapping subhaloes are 
still rare, and because it is not clear whether these pairs have 
mixed completely, we do not integrate them together in the 
current version of HBT. We also notice that one SUBFIND 
subhalo with 2 x 10"^ particles is missed by HBT, indicating 
the tracing algorithm is still not perfect although it is al- 
ready massively optimized. However, as statistically found 
in our comparison using the cosmological simulation, this 
kind of missed subhaloes can be as large as 10** particles for 
SUBFIND, amounting to as high as 10 percent in number 
near that mass, while HBT only misses less than 1 percent 
for the highest missed mass of 10^ particles. 

A second situation is when one subhalo finder only iden- 
tifies the most inner part of a subhalo and assigns a majority 
of surrounding bound particles to the central subhalo, while 
its correspondence identifies both parts. When doing the 
match, the surrounding particles from the correspondence 
outweighs the inner part, linking the correspondence to the 
central subhalo. This situation can be found as two adjacent 
perpendicular lines connecting two unilateral subhaloes and 
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the central. It happens for both finders. In the SUBFIND 
case, it is due to the truncation around the radius of equal- 
ity as illustrated in section 13.41 while in the HBT case we 
attribute this to the effect of local accretion discussed in 
Appendix [B] 

The third situation is similar to the second one but 
happens for binary subhaloes, which has been discussed in 
the main text as "sub-in-sub" switch. 



APPENDIX E: ACCURACY OF TWO 
ESTIMATORS FOR THE SUBHALO MASS 
FUNCTION 

The subhalo mass function /(m) = dN/dlnm where m is 
subhalo mass (either Msub or Msub/Mhost) can be estimated 
directly from two discrete estimators as 

AiV 



fa[m) = 



fb{m) 
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A Inm 

mi <m^ <m2 
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where mi and m2 are the lower and upper limits of subhalo 
mass bin, AA^ is the number of subhaloes in the bin and m 
is the average subhalo mass in the bin. In the limit Am 0, 
we have m = m and /a(m) — fbijn) = f(m). For a finite 
mass bin, suppose /(m) = ignoring the normalization, 
then 
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Figure El. Bias of mass function estimators as a function of 
mass bin size Alog(m) = log(fc). 
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It is ready to see that the estimator fa recovers exactly /. 
For a more general form /(m) — m let k = m2/mi, then 
the biases of the two estimators are given as 
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Since the biases do not depend on m, the slope of the mass 
function is recovered for both estimators while the normal- 
izations differ. In Figure [eTI we show the bias of these two 
estimators for 7 — 0.9. With bin size Alog(m) — 0.5, ft will 
underestimate the mass function by 10 percent. Throughout 
this work we use fa to calculate subhalo mass function ex- 
cept in Figure [3] and Figure [5] where fb is used to give the 
total subhalo mass contained in each mass bin. 
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